
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Test: Constancy of the Slope Coefficients      %
%                Monte Carlo Experiment                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

L2=100; % number of auctions with N=2
L3=100; % number of auctions with N=3
L=L2+L3;
N=[2,3];
alpha=0.50:0.10:0.80; % quantile levels alpha
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];

nsim=10000; % Number of simulations
B=500; % Bootstrap replications
toler=0.1; % Parameter of tolerance for estimation
m=size(alpha,2);
Rej1=0; % auxiliary sum for rejection probability
Rej5=0; 
Rej10=0;

% Simulation 

for i=1:nsim
        
        % Data Generating Process
        eps2=exprnd(1,L2*N(1),1); % Regression Errors N=2
        eps2=reshape(eps2,N(1),L2);
        eps2=sort(eps2,'descend');
        eps2=eps2(2,:);
        eps3=exprnd(1,L3*N(2),1); % Regression Errors N=3
        eps3=reshape(eps3,N(2),L3);
        eps3=sort(eps3,'descend');
        eps3=eps3(2,:);
        
        X2=normrnd(0,1,L2,1); % Covariates N=2
        X3=normrnd(0,1,L3,1); % Covariates N=3
        x2=[ones(L2,1),X2];
        x3=[ones(L3,1),X3];
        x=[x2;x3];
        p=size(x,2); 
        
        s2=1+(0*X2);
        s3=1+(0*X3);
        
        gamma0=zeros(1,m); % Initial guess for the intercept
        gammap=1; % Initial guess for the slope
        gamma = [0,1]';
        
        wb2=0+gammap*X2+s2.*eps2'; % Winning bids N=2
        wb3=0+gammap*X3+s3.*eps3'; % Winning bids N=2
        wb=[wb2;wb3];
        
        % Computation of the Test Statistic for Constancy of Slope Coefficient    
        [obj_cqr,gamma0,gammap,res_cqr]=CQR(psi, p, wb, x, toler,gamma0,gammap); % Estimation under H0
        obj0=sum(obj_cqr)/m;

        for k=1:m
            [obj,res,~]=MM(psi(:,k), p, wb, x, toler, gamma); % Estimation under H1
            resk(:,k)=res;
            objk(k,:)=obj;
        end
        obj1=sum(objk)/m;
        M=obj0-obj1; % test statistic
    
        % Random Weighting Bootstrap of the Test Statistic for Constancy of Slope Coefficient 
        
    for b=1:B
        
        w = [mnrnd(L2,ones(L2,1)/L2)';mnrnd(L3,ones(L3,1)/L3)']; % weights
         
        [obj_cqr_wb,~,~,~]=CQR_w(psi, p, wb, x, w, toler,gamma0,gammap); % Estimation under H0
        obj0wb=sum(obj_cqr_wb)/m;
        
        for k=1:m     
            res_cqrw=res_cqr(:,k).*w; % residuals under H0 weighted by w
            resw=resk(:,k).*w; % residuals under H1 weighted by w    
            obj_cqrw(k,:)=sum(psi(:,k).*res_cqrw)-sum(res_cqrw(res_cqrw<0)); % objective function under H0 weighted by w
            objw(k,:)=sum(psi(:,k).*resw)-sum(resw(resw<0)); % objective function under H1 weighted by w
            % Weighted objective functions
            [obj_wb,~,~]=MMw(psi(:,k), p, wb, x, w, toler, gamma);
            objk_wb(k,:)=obj_wb;
        end
        obj0w=sum(obj_cqrw)/m; 
        obj1w=sum(objw)/m;
        Mw=obj0w-obj1w; % test statistic
        
        obj1_wb=sum(objk_wb)/m;
        Ms=obj0wb-obj1_wb;
        Mstar(:,b)=Ms-Mw; % bootstraped test statistic
        
    end
    
    % Critical Values
    c99=quantile(Mstar',0.99);
    c95=quantile(Mstar',0.95);
    c90=quantile(Mstar',0.90);
    
    % Computation of rejection probability 
    if (M>c99)
    Rej1=Rej1+1;
    end

    if (M>c95) 
    Rej5=Rej5+1;
    end

    if (M>c90)
    Rej10=Rej10+1;
    end
    rej=Rej5/nsim;
    i
end

% Rejection probability
Rej1=Rej1/nsim
Rej5=Rej5/nsim
Rej10=Rej10/nsim



